
#include <cassert>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;

constexpr double vel{34.02970287557621}, mu{1.7892976260350732e-05},
    ref_density{1.22498}, cspeed{340.297};

struct Fun {
  double value = 0.0;
  double operator()(double x) {
    double tmp = (37.0 / 315 - x / 945 - x * x / 9027);
    return tmp * tmp * x - value;
  }
};

struct FK {
  double operator()(double x) {
    return 2.0 * (37.0 / 315 - x / 945 - x * x / 9027) *
           (2.0 - 116.0 / 315 * x + (2.0 / 945 + 1.0 / 120) * x * x +
            2.0 / 9027 * x * x * x);
  }
};

template <typename T> struct Bisection {
  double operator()(T &f, double a, double b, const double tol = 1.0e-6) {
    double fa, fb, fc;
    fa = f(a);
    fb = f(b);
    double c = 0.5 * (a + b);
    fc = f(c);

    while ((fabs(fc - 0) > 1.0e-15) && (c - a > tol)) {
      if (fa * fc < 0.0)
        b = c;
      else
        a = c;

      c = 0.5 * (a + b);
      fc = f(c);
    }

    return c;
  }
};

int main() {

  vector<double> x, y, u, v, density_vec;
  ifstream fdata{"surface_flow.dat"};
  int num_points;
  fdata >> num_points;
  x.reserve(num_points);
  y.reserve(num_points);
  density_vec.reserve(num_points);

  u.reserve(num_points);
  double u1, u2, tmp_d;
  for (int i = 0; i < num_points; ++i) {
    fdata >> x[i] >> y[i] >> density_vec[i] >> u1 >> u2;
    // 将无量纲量转换为有量纲量
    u[i] = sqrt(u1 * u1 + u2 * u2) * cspeed / density_vec[i];
    density_vec[i] *= ref_density;
    // cout << u[i] << endl;
  }
  fdata.close();
  // 计算占位  没有归一化
  vector<double> rad_len(num_points);
  rad_len[0] = 0.0;
  double tmp{0.0};
  for (int i = 0; i < num_points - 1; ++i) {
    tmp += sqrt((x[i + 1] - x[i]) * (x[i + 1] - x[i]) +
                (y[i + 1] - y[i]) * (y[i + 1] - y[i]));
    rad_len[i + 1] = tmp;
    // cout << rad_len[i] << "\n";
  }

  // 计算梯度  将不均匀看做均匀处理  网格需均匀过渡
  double u0_grd, u0_ggrd;
  u0_grd = (4 * u[1] - 3 * u[0] - u[2]) / (rad_len[2] - rad_len[0]);
  // u0_grd = (u[1] - u[0]) / rad_len[1];
  u0_ggrd = 4 * (u[0] - 2 * u[1] + u[2]) /
            ((rad_len[2] - rad_len[0]) * (rad_len[2] - rad_len[0]));
  // u0_ggrd = 0.0;
  // u0_grd = 0.0;

  cout << "u0_grd = " << u0_grd << "\n"
       << "u0_ggrd = " << u0_ggrd << endl;

  double z0, z0_grd;
  z0 = 0.0770 / u0_grd;
  z0_grd = -0.0652 * u0_ggrd / (u0_grd * u0_grd);

  // 积分
  vector<double> delt1(num_points, 0.0), delt2(num_points, 0.0),
      delt(num_points, 0.0), tau0(num_points, 0.0);

  Fun K;
  FK fk;
  Bisection<Fun> bisec;
  double z1, z2, z_grd, u_grd, k1, gam, alpha2;
  z1 = z0;
  z_grd = z0_grd;

  int i;
  vector<double> gamma(num_points, 0.0);
  for (i = 1; i < num_points - 1; ++i) {
    z1 = z1 + z_grd * (rad_len[i] - rad_len[i - 1]);

    u_grd = (u[i + 1] - u[i - 1]) / (rad_len[i + 1] - rad_len[i - 1]);
    // u_grd = (u[i + 1] - u[i]) / (rad_len[i + 1] - rad_len[i]);
    K.value = z1 * u_grd;
    cout << "k_value = " << K.value << endl;
    gam = bisec(K, -12.00001, 12.00001);

    if (gam > 12.0 || gam < -12.0)
      break;
    gamma[i] = gam;
    cout << "gam = " << gam << endl;
    delt2[i] = sqrt(z1 * mu / density_vec[i]);

    alpha2 = 1.0 / 63 * (37.0 / 5 - gam / 15 - gam * gam / 144);
    delt[i] = delt2[i] / alpha2;
    delt1[i] = (0.3 - gam / 120) * delt[i];
    tau0[i] = (2.0 + gam / 6) * mu * u[i] / delt[i];

    z_grd = fk(gam) / u[i];
  }
  gam = 7.052;
  delt2[0] = sqrt(z0 * mu / density_vec[0]);
  alpha2 = 1.0 / 63 * (37.0 / 5 - gam / 15 - gam * gam / 144);
  delt[0] = delt2[0] / alpha2;
  delt1[0] = (0.3 - gam / 120) * delt[0];
  tau0[0] = (2.0 + gam / 6) * mu * u[0] / delt[0];

  // cout << "i = " << i - 1 << endl;
  //
  double length = 0.0;
  for (int i = 0; i < num_points - 1; ++i) {

    length += sqrt((x[i + 1] - x[i]) * (x[i + 1] - x[i]) +
                   (y[i + 1] - y[i]) * (y[i + 1] - y[i]));
  }

  ofstream fout{"solution.dat"};
  fout << scientific << setprecision(15);
  fout << "VARIABLES = \"length\", \"delta\", \"delta1\", \"delta2\", "
          "\"tau0\"\n";
  fout << "ZONE  I=" << i - 1 << endl;
  // 输出计算结果
  for (int j = 0; j < i; ++j) {
    fout << rad_len[j] / rad_len[num_points - 1] << "  " << delt[j] << "  "
         << delt1[j] << "  " << delt2[j] << "  " << tau0[j] << "\n";
  }
  fout.close();

  fout.open("gam_re.dat");
  fout << scientific << setprecision(15);
  fout << "VARIABLES = \"gam\", \"delta\"\n";
  fout << "ZONE  I=" << i - 2 << endl;

  for (int j = 1; j < i; ++j) {
    fout << gamma[j] << "  " << u[j] * delt1[j] * density_vec[j] / (mu) << "\n";
  }
  fout.close();

  // test
  //   Fun K;
  //   FK fk;
  //   Bisection<Fun> bisec;
  //   K.value = -0.1474;
  //   double tmp_gam = bisec(K, -12.0001, 12.0001);
  //   cout << tmp_gam << "\n" << fk(tmp_gam) << endl;

  //   K.value = 0.0;
  //   double tmp_x;
  //   for (int i = 0; i < 1000; ++i) {
  //     tmp_x = (i - 500) * 0.03;
  //     cout << tmp_x << "  " << K(tmp_x) << "\n";
  //   }

  return 0;
}
